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The dune skeleton model is a reduced model to describe the formation process and dynamics of 
characteristic types of dunes emerging under unidirectional steady wind. Using this model, we study 
the dependency of the morphodynamics of transverse dunes on the initial random perturbations and 
the lateral field size. It was found that i) an increase of the lateral field size destabilizes the transverse 
dune to cause deformation of a barchan, ii) the initial random perturbations decay with time by the 
power function until a certain time; thereafter, the dune shapes change into three phases according 
to the amount of sand and sand diffusion coefficient, iii) the duration time, until the transverse dune 
is broken, increases exponentially with increasing the amount of sand and sand diffusion coefficient. 
Moreover, under the condition without the sand supply from windward ground, the destabilization 
of transverse dune in this model qualitatively corresponds to the subaqueous dunes in water tank 
experiments. 
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1. INTRODUCTION 

Erosion due to wind sculpts deserts on Earth and sur- 
faces on Mars, Titan into sand dunes such as barchan, 
transverse, longitudinal, star-shaped, and dome-shaped 
[ll-Q. As the dominant factors dictating several dune 
shapes, the steadiness of wind direction and the amount 
of available sand in each dune field are known For 
example, a unidirectional wind generates barchans of a 
crescent-shaped or transverse dunes extending perpen- 
dicular to the wind direction depending on the amount 
of available sand, whereas a bidirectional wind generates 
longitudinal dunes extending parallel to the sum of two 
wind directional vectors. Additionally, a multidirectional 
wind generates star dunes consisting of multiple crest 
lines extending from the top to several directions. 

As one of remarkable subject in dune studies, the sta- 
bility of transverse dunes has remained to be an open is- 
sue for both geomorphology and geophysics. In rescaled 
water tank experiments, an isolated transverse dune was 
shown unstable and deformed into a set of barchans as 
the end-stage Also, computer models have repro- 

duced the qualitative and quantitative morphodynamics 
similar to subaqueous dunes [IT14T71 ]. Recently, theoreti- 
cal approaches for the stability of transverse dunes have 
been conducted by the methodology considered the mass 
conservation of sand and the sand flows on dunes [18l421| . 
Especially, Niiya et al. consider the dynamics of a dune 
as a combination of two-dimensional cross sections (here- 
after, 2D-CSs) parallel to the wind direction and pro- 
posed the dune skeleton model (hereafter, DS model) con- 
sisting of coupled ordinary differential equations. This 
model is based on several assumptions used in the pre- 
vious analytical model named "aeolian/aqueous barchan 
collision dynamical equations (ABCDE)" by Katsuki and 
Nishimori. ABCDE is able to describe the collision dy- 
namics of two 3D barchans focusing on the central 2D-CS 
of barchans [22j, |23[ . 

So far, the DS model has successfully reproduced three 
typical shapes of dunes, straight transverse dune, wavy 



transverse dune, and barchan, depen ding on the amount 
of available sand and wind strength [l8|. Moreover, the 
reduced DS model, which is a further simplified DS model 
with a two-variable ordinary differential equation, en- 
abled the elucidation of the mechanism of transition be- 
tween different dune shapes using a bifurcation analysis 
[2ll ] . However, in these previous studies, the initial con- 
dition was given as a sinusoidal curve with small ampli- 
tude and single-wavelength. Thus, the stability of dunes 
against random perturbations is yet to be identified. In 
this study, we investigate the effect of initial perturba- 
tions on the stability of transverse dunes, and the field 
size effect on the stability is also studied. 



2. DUNE SKELETON MODEL 

The DS model covers the formation processes of 
barchans and transverse dunes, both of which are gen- 
erated under a unidirectional steady wind. This model is 
roughly based on three considerations. First, the dunes 
consist of 2D-CSs as mentioned in the previous section. 
Second, a lateral distance between neighboring 2D-CSs 
is set constant. Third, a combination of two forms of 
sand movement, intra-sand movement within each 2D-CS 
and inter-sand movement between neighboring 2D-CSs, 
is considered to govern the macroscopic morphodynamics 
of dunes. 

3D barchans are isolated on hard ground in both wind 
and lateral directions, whereas 3D transverse dunes are 
not isolated because they extend in the lateral direc- 
tion. However, in the wind direction, the laterally ex- 
tending 3D transverse dunes are assumed to be sepa- 
rated by inter-dune hard ground. In addition, consid- 
ering the observation that the 2D-CSs of barchans and 
transverse dunes very roughly show a scale-invariant tri- 
angular shape, we assume that the angles of their upwind 
and downwind slopes (9 and (p, respectively as shown in 
Fig. [l} are maintained constant through their migration, 
irrespective of their size. From these assumptions, the 
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FIG. 1. Intra-2D-CS sand flow. The over-crest sand flux q 
and sand trapping efficiency Te govern the intra-2D-CS flow. 
The green solid and dashed lines indicate sand deposition on 
the downwind slope and escaping sand to the leeward ground, 
respectively. 
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FIG. 2. Inter-2D-CS sand flow J u (upwind) and 3d (down- 
wind). Both the red and blue arrows indicate local sand 
transport depending on the height difference. The red and 
blue areas correspond to J u and Jd, respectively, by summing 
the local sand transport from the slope's foot to the 2D-CS's 
crest. 



horizontal (i.e., wind directional) position and the height 
of each 2D-CS are uniquely determined if only the coor- 
dinate (x, h) of its crest is given. Moreover, the empirical 
geometrical constants A, B, and C are introduced as 



tan 9 tan w 
A = — ,B 



tan ip 



C 



tan 9 



tan 9 + tan ip tan 9 + tan ip tan 9 + tan tp 

where A, B, and C are set to 1/10, 4/5, and 1/5, re- 
spectively, reflecting the typical 2D-CS profiles of real 
barchans and transverse dunes. 

As mentioned above, sand flow is classified into two 
forms: (a) the intra-2D-CS flow and (b) the inter-2D-CS 
flow. The intra-2D-CS flow along the upwind slope is 
uniquely determined if the over-crest sand flux q and the 
incoming sand flux from the windward ground f m are 
given (Fig. [T]). The over-crest sand flux q determines the 
erosion rate of the 2D-CS's upwind slope and the over- 
crest blown sand deposits along the downwind slope or 
directly escapes to the leeward inter-dune ground. The 
deposition ratio T E of the over-crest sand along the down- 
wind slope is assumed as an increasing function of height 
of the specific form 



T E (h) 



h 



1.0 + h 



(1) 



Note that T E (h) is termed as the sand trapping efficiency 
[2j| and Eq. ([T]) roughly represents the case of typical 
shear velocity w* = 0.4 m s _1 at the dune crest, where h is 
the height in a meter unit. The quantity of q reflects the 
over-dune wind strength; here, we assume q, where 0.1 < 
q < 1-0, to be constant in each desert field, independent 
of the 2D-CS's height. In addition, all of the incoming 
sand flux from the windward ground, f m , deposits on 
the upwind slope of dunes. 

The inter-2D-CS flow J u u_^j\ / J^fj^i) occurs only be- 
tween the upwind/downwind slopes of the neighboring 
2D-CSs, i and j (Fig. [2]). Locally, most of the lateral 
sand transport is determined by the height difference be- 
tween neighboring 2D-CS's slopes. Therefore, we assume 
that the total flux is the sum of local sand transport 
from the slope's foot to the 2D-CS's crest. Namely, the 



flux is roughly considered as the lateral diffusion depend- 
ing on the height difference, though the consideration of 
the overlap length of slopes causes a nonlinearity in the 
present form of the inter-2D-CS flux. The specific forms 
of J u (i^j) and Jd{j~>i) are 
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where Aw, set as Aw = 1 hereafter, is the lateral inter- 
val between neighboring 2D-CSs. These quantities cor- 
respond to the colored areas in Fig. [2] multiplied by dif- 
fusion coefficients. Here, the upwind and downwind dif- 
fusion coefficients (D u and Dd, respectively) control the 
amount of inter-2D-CS sand flow on the respective sides 
of the slopes and reflect the over-dune wind strength. 

With consideration of the above intra- and inter- 
sand flows, the dynamics of the coordinates (xi,hi) (i = 
1, • • • ,N) of the 2D-CS's crest are given as a system of 
coupled ordinary differential equations [H, : 



dxi 1 
dt hi 



dhi _ A 
dt hi 



q(BT E (hi) + C) + (B-fyj-rt) + CJ*.(i_tf)) - °/i 
i=i±l 

(3a) 



q(TE(hi) — 1) + (J d (j~>i) — Ju(i^j)) + fl 
j=i±i 



(3b) 



We also introduce the annihilation rule of the 2D-CSs at 
the lateral edges of a dune; this rule is required to simu- 
late the shrinking process of dunes. This rule is applied 
in the cases where hi is lower than by the escaping of 
sand from i-th 2D-CS (Fig. [3][a)) or where the overlap 
between the 2D-CS at the edge and its neighbor vanishes 
(Fig. [3Jb)). If the annihilation rule is applied to i-th 
2D-CS, we consider that the state of i-th lane gets cor- 
responding to the exposed state of ground, that is, i-th 
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FIG. 3. Schematic of annihilation rule, (a) Annihilation by 
escaping of sand from 2D-CS. h = mean the exposed state 
of ground, (b) Annihilation by leaving 2D-CS from dune. The 
blue area is overlap between the 2D-CS at the edge and its 
neighbor. 
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FIG. 4. Typical dune shapes depending on Gap(t) by nu- 
merical simulation: (a) straight transverse dune, (b) wavy 
transverse dune, and (c) barchan. 



2D-CS is taken out of the calculation though the virtual 
crest of height zero is fixed at the foot of downwind slope 

as 



hi = 0, 



Xi = max 
j=i±i 



In order to compare the results obtained from DS 
model with the real desert dunes, we estimate the value 
of q, D u , and in physical units on the basis of the 
observed data in real dune fields. The migration velocity 
of single 2D-CS in our model is assumed as 



dx q . . 

H = h TE{h) = T 



(I 



h 



(4) 



Indeed, in field observations [25-27], the migration speed 
c (m/year) and size H (m) of barchans are reasonably 



fitted as [28( 



Q Q/Hp 

H + H 1 + H/H ' 



(5) 



where Q is the volumic sand flux (m 2 /year), Ho is the 
cut-off height ranging from 1.8 m (Finkel 1955-58) to 10.9 
m (Hastenrath 1958-64). Comparing Eq. (gj) with Eq. 
([5]), the relationship between our model and real data is 



given as 



h=— q=^- = c 
Ho Ho 



1 + # 

-no 



Here, we roughly assume that a dune with 5 m height 
migrates 10 m per year and the cut-off height is 1 m, 
namely, c = 10, H = 5, and Ho = 1 in meter-year unit. 
Converting time unit from year to day, we get about q « 
0.164 m 2 /day. Also, q used in our model takes the value 
from 0.1 to 1.0. Therefore, characteristic time unit and 
space unit accompanied with physical values used in our 
model are approximately estimated as a day and a meter, 
respectively. 



RESULTS 



Numerical simulation of Eqs. (|3aj) and (|3b[) is carried 
out using several numbers of 2D-CSs, N — 200, 600, 



1000, 1400, and 1800. The lateral boundary condition 
is set as periodic, whereas the wind directional boundary 
condition is set such that sand escaping from the leeward 
boundary is redistributed uniformly from the windward 
boundary. Thus, mass conservation, i.e., the total area of 
2D-CSs, is maintained throughout the simulation unless 
the annihilation of a 2D-CS occurs. Moreover, according 
to the rule for the redistribution of escaping sand, the 
incoming sand flux to each 2D-CS from the windward 
ground f m is expressed as 



1 ^ , 



T E {hi)), 



i=l 



where Te is the sand trapping efficiency given by Eq. 
([I]). In order to elucidate the stability of the transverse 
dunes against a perturbation, the initial dune shape is set 
slightly fluctuated from a straight transverse dune. More 
specifically, while the initial heights of 2D-CS's crests are 
set uniform, i.e., hi(0) = Ho, the fluctuation of initial 
wind directional positions of the crest of each 2D-CS is 
set random with an amplitude depending on Ho', that 
is, £Ej(Q) = e(Ho). Here, e(H ) is a random number uni- 
formly distributed between and Hq/IO. 

Like in the previous paper (Niiya et al., 2010), we clas- 
sify the obtained dune shapes into three different phases 
according to the quantity, 



Gap(t) 



1 



N 



(6) 



Gap describes the degree of sinuosity for a transverse 
dune (Fig. gj. 

I) : If Gap(l0 8 )/Gap(0) < 10~ 2 is satisfied, we regard a 

straight transverse dune is formed. The phase in 
this state is referred to as the "ST phase" . 

II) : If Gap(10 s )/Gap(0) > 10~ 2 is satisfied, we regard a 

wavy transverse dune is formed. The phase in this 
state is referred to as the "WT phase" . 

III) : If the annihilation of a 2D-CS occurs before t = 

10 8 , we regard the wavy transverse dune deforms 
to the barchan. The phase in this state is referred 
to as the "B phase" . 
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FIG. 5. Phase diagram of dune shapes obtained from numer- 
ical simulations by varying the two pairs of control parame- 
ters: (a) (Ho,N) and (b) (D U ,N). The blue, green, and red 
marks denote the ST, WT, and B phase independent of the 
initial random perturbation, respectively. Other parameters 
are set as (a) D u = 0.1 and (b) Ho = 30. (c) Typical shape of 
transverse dune against field size: (red) N = 1000 and (blue) 
TV = 1800. The dune with red crest line represents the sta- 
tionary wavy transverse dune in WT phase, whereas the dune 
with blue one represents the dune at the moment when the 
annihilation of a 2D-CS occurs. The parameters are fixed as 
(#o,A0=(20,0.1). 



Note that, through an analytical approach of reduced DS 
model with only two 2D-CSs, the transverse dunes are 
shown to be attracted to one of ST, WT, and B phases 

m. 

In this simulation of the DS model, we fix the over- 
crest sand flux in each 2D-CS, q — 0.5 (m 2 /day), which 
reflects the wind strength in a dune field, downwind diffu- 
sion constant given in Eq. (|2b|) is set Dd = 0.1 (m 2 /day), 
and the lateral period between neighboring 2D-CSs is as- 
sumed as 1 m. Then, three parameters, N, Hq, and D u 
are varied, where N (m) is the number of 2D-CSs, i.e., 
the lateral field size in meter, Hq (m) is the uniform ini- 
tial height in meter of 2D-CSs, which controls the total 
amount of sand in the system, and D u (m 2 /day) is the 
upwind diffusion coefficient introduced in Eq. (|2a[) . The 
numerical calculations are performed until the simulation 
time reaches t = 10 s (day) or the annihilation of a 2D-CS 
occurs. If annihilation occurs, a dune rapidly shrinks and 
eventually disappears. We take five ensembles of calcula- 
tion with different amplitude of initial perturbations for 
each set of parameters. Obtained results are as follows. 

First, according to the above classification rules 
I)~III), we verify the influence of the lateral field size 
on the final dune phase. In Figs. EJa) andOJb), the fi- 
nal dune phases obtained in two sets of parameter spaces 
(H ,N),(D U ,N) are classified. Figure [5jc) shows the 
lateral system size dependency of the stability of trans- 
verse dunes, where initially same height (Hq = 20) of 
transverse dunes are set in different lateral size (N — 



FIG. 6. Typical time evolutions of dune sinuosity using Gap: 
(a) Hq and (b) D u . The solid lines are selected from the simu- 
lations with N — 1000 in Fig. [5ja). a is the power index, (c) 
Comparison DS model and one-dimensional diffusion equation 
until t = 10 . (d) Comparison Gap (solid) and height differ- 
ence h (dashed) in the system. The lines are picked from (a): 
(green) Ho = 19 and (purple) Ho = 30. 



1000, 1800) of fields. The transverse dune in a longer lat- 
eral size of field will break, whereas the transverse dune 
in a shorter lateral size of field is kept steadily in the 
WT phase. It means that the critical wave length A c for 
the occurrence of the instability of transverse dune exists 
between 1000 m and 1800 m in this case. It should be 
noted that A c depends also on D u and Dd, so it is required 
to estimate the specific values of them in respective dune 
fields, which is the next issue. Another relevant factor for 
the stability of a transverse dune in the present model is 
the effect of homogeneous redistribution of sand from the 
upper boundary, that seems to stabilize a transverse dune 
through supplying a same amount of sand at the upwind 
slopes of both the thin-and-low part and the thick-and- 
high part of a transverse dune, the slowing down effect 
for the former part is more obvious than that for the 
latter part. Therefore, the destabilization by the fasten 
migration of the lower part is suppressed. 

Second, in order to see the evolution process of trans- 
verse dunes with a fixed lateral size (N — 1000), we fo- 
cus on the time evolution of Gap. Specifically, we pick 
up four points from Fig. O^a): (H ,N) = (11,1000) in 
B phase (named point-a), (H ,N) = (19,1000) in WT 
phase near the boundary to B phase (named point-b), 
(H , N) = (21, 1000) in WT phase near the boundary to 
ST phase (named point-c), and (H ,N) = (30,1000) in 
ST phase (named point-d). Figures EJa) andEtb) show 
the time evolution of Gap, and it is characteristic that 
the initial fluctuation of all above defined points decays 
with time in the power of a ~ —1.5 until t = 10 5 ; that 
is, the initial transverse dunes temporarily approach the 
straight shape of transverse dune. Thereafter, Gap at 
point-a turns into a drastic increase, where the annihi- 
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lation of a 2D-CS takes place, namely, the breaking of a 
transverse dune occurs. In contrast, Gap at the remain- 
ing three points, b, c and d, once increases the decreasing 
rate, among which Gap at point-b and point-c turn into 
increasing after certain durations of time, and finally con- 
verge to the respective final states. Especially, at point-b, 
Gap oscillates with a small amplitude in the final state, 
whereas Gap settles to a constant value at point-c. 

To understand these characteristic evolution of Gap, 
we divide the simulation period into two parts; i) power- 
law decaying time region until t — 10 5 and ii) drastic 
growing or decaying time region after t = 10 5 . For i), 
remember that the perturbation 5xi added to the ini- 
tial straight transverse dune with x% = xq , hi = Hq was 
randomly selected value between and Hq/10. Here, 
we assume that the power decay of Gap corresponds to 
the smoothing process of the initial perturbation caused 
by the lateral sand transport. To check this, the one- 
dimensional diffusion equation 



dfl - r>( f 
~dt ~ D{fl+1 



/i-l-2/0 (» = 1,2,..-,L) (7) 



under the same periodic condition as the DS model is 
numerically calculated, where the variables in Eq. ([7]) 
correspond as (fi,L,D) « (xi,N,D u ) with /i(0) = 
f(0)+6fi (Sfi = Sxi), and the time evolution of Gap(t) = 

Si^i (/»(*) — /i+i(0) 2 i s measured. Then, the decay rate 
of Gap in the discrete diffusion equation (Eq. 0) is con- 
firmed to obey the power-law with exponent a ~ —1.5. 
This exponent is analytically derived from a simple di- 
mension analysis of the diffusion equation (Fig. Htc)). In 
this way, the power-law decay of Gap in the DS model in 
the time region i) is strongly guessed to originate from 
the lateral diffusion of sand. 

For ii), the time evolution of Gap after t = 10 5 shows 
the exponential growth or decay (Fig. EJd)), where the 
amplitude of single wave shown in Fig. [SJc) increases or 
decreases. In order to see this drastic change of the time 
evolution of Gap from another aspect, we measure the 
difference between the maximum and minimum heights 
of 2D-CSs, h = h max — h m i n . Here, Gap increases fol- 
lowing the increase of h because the velocity of 2D-CSs 
is inversely proportional to their height from Eq. (|3a|) . 
Actually, in Fig. EJd), h turns into increasing before the 
increase in Gap. Thus, the transition of Gap from decay 
to growth is induced by the increase in height difference 
of 2D-CSs. In the next, we consider the mechanism of the 
exponential increase or decrease in Gap. In the diffusion 
equation (Eq. ([?[))■ if the wave number k of transverse 
dune is determined, f%(t) decay exponentially with time 
according to the decay rate —Dk 2 . In contrast, the DS 
model consists of the nonlinear coupling of position and 
height of 2D-CSs as seen in Eqs. (|5a|) and (I3bp . which 
may act as the factor to remain the wavy transverse dune 
stable. Namely, the exponential decay of Gap is gener- 
ated by the diffusion, whereas the exponential growth 
comes from the coupling between position and height. 
This problem is left as the next issue. 
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FIG. 7. Duration times of transverse dunes and growth rate 
of Gap after diffusion process: (a, c) Ho and (b, d) D u . (3 
is the exponent, (a, b) The colors of the points indicate the 
different N used for the simulations. Each point is the average 
duration time per simulation, (c, d) The fitting values by 
exponential function. TV is fixed as TV = 1800. 



Finally, we study the duration time T d , which rep- 
resents the time interval of transverse dunes before the 
break of the crestline. Note that the maximum value of 
Td equals the simulation time t = 10 s . Figures [TJa) and 
[Tfb) show that the values Td lower than 10 8 obtained in 
the B phase is nicely fitted by exponential function of two 
control parameters as Td oc exp(0.2i?o + 4SD U ), though 
an increase in TV extends the range of the B phase. To 
elucidate the cause of these exponential increases of Td 
as functions of Hq and D u , we focus on the exponen- 
tial growth region of Gap after diffusion process in Fig. 
[6l where Gap{t) oc exp(crt) holds with growth rate a. 
If Gap exceeds a certain value and the transverse dune 
breaks, therefore, the duration time is roughly estimated 
as Td oc 1/er. We measure a from simulation data as 
shown in Figs. [Tfc) and [7Jd) and a in the B phase is 
nicely fitted as a oc exp(— Q.2Hq — 48D U ), which rela- 
tions correspond to the dependency of T<j on Hq and D u 
shown in Figs. [TJa) and[7]Jb). 



4. CONCLUSION 

In this study, we discussed on the stability of transverse 
dunes against random perturbations and its dependency 
on lateral size using the DS model. Numerical simula- 
tions show that i) an increase in the number of 2D-CSs, 
namely, the lateral size increase, destabilizes the trans- 
verse dunes, ii) in the initial stage of simulations, the 
small scale of fluctuation of a transverse dune is smoothed 
by the lateral sand diffusion, and iii) after the initial dif- 
fusion regime of time, transverse dunes evolve according 
to the specific growth rate depending on two control pa- 
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FIG. 8. Deformation of transverse dune without redistribu- 
tion rule of sand. 

rameters, Hq and D u . The uniform redistribution rule 
of sand in the present model affects on the stability of 
transverse dunes. Specifically, the escaping sand from 
a transverse dunes is uniformly supplied from the wind- 
ward boundary in the present model. 

Note that most of previous water tank experiments and 
mathematical models have indicated that the inevitable 
instability of transverse dunes causes the deformation of 
them into barchans. In such cases, the boundary condi- 



tion in wind direction was set open, namely, no-sand was 
supplied from the windward ground. This boundary con- 
dition is definitely different from that in the present DS 
model. If the open boundary is used in the present DS 
model, the total amount of sand consistently decreases 
and the transverse dune can not keep the stable states. 
To demonstrate this, we conduct the numerical simula- 
tion without the redistribution rule of sand under the 
same initial condition as the above mentioned case. The 
initial condition is randomly given as with this paper. 
Then, the initially straight transverse dune deforms into 
the wavy shaped dune such like barchanoid, and further 
simulation causes the emergence of many barchans, even 
though the barchans eventually disappear (Fig. |8|). In 
this way, under the no-sand supply condition, our model 
shows the instability of transverse dunes similar to pre- 
vious experiments and models. Based on the this re- 
sult, it is reasonably guessed that the morphodynamics 
of dunes vary largely depending on the supplying condi- 
tion of sand. 
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